sb urban tree parameters

Five different tree species were used to select parameter sets to represent different types found around Santa Barbara.

1. coastal live oak Quercus agrifolia - native broadleaf evergreen
2. canary island pine Pinus canariensis - needleleaf evergreen
3. victorian box Pittosporum undulatum - nonnative broadleaf evergreen
4. blue gum eucalyptus Eucalyptus globulus - nonnative broadleaf deciduous
5. california sycamore Platanus racemos - native broadleaf deciduous

Steps taken to select the parameters were as follows:

Plots of model output

## summarized annually 
dfyr %>%  
  ggplot() + geom_line(aes(x=as.factor(year), y=lai, col=code, group=code))

dfyr %>% dplyr::filter(year>2009 & year < 2020) %>% 
  ggplot() + geom_line(aes(x=as.factor(year), y=lai, col=code, group=code))

dfyr %>% 
  ggplot() + geom_line(aes(x=as.factor(year), y=plantc, col=code, group=code))

#dfday %>% ggplot() + geom_line(aes(x=yd, y=lai, col=code))

dfyr2 %>% ungroup() %>%
  dplyr::filter(year ==2011 | year==2014 | year==2017) %>% 
  ggplot() + geom_boxplot(aes(x=code, y=lai, fill=as.factor(year))) + ggtitle("annual avg LAI") 

dfyr2 %>% ungroup() %>%
  dplyr::filter(year ==2011 | year==2014 | year==2017) %>% 
  ggplot() + geom_boxplot(aes(x=code, y=plantc, fill=as.factor(year))) + ggtitle("annual avg plantc")

Plots of rs data and model outputs

# load in mikes data 
rs_dat <- read.csv("../../data/alonzo_data_df.csv")
rs_dat <- dplyr::filter(rs_dat, fractional_cov > 0.7)
trees_rs <- rs_dat %>% dplyr::filter(
    class==81 | #QUAG 
    class==74 | #PLRA
    class==73 | #PIUN
    class==66 | #PICA
    class==34 ) #EUGL
## remove outliers 
trees_rs_rm <- trees_rs[!trees_rs %in% boxplot.stats(trees_rs$carbon)$out,]

ggplot(trees_rs_rm) + geom_boxplot(aes(x=as.factor(class), y=lai, fill=as.factor(class))) + ggtitle("LAI estimates from alonzo 2016")

## Warning: Removed 2143 rows containing non-finite values (stat_boxplot).

trees_rs_rm$code <- NA
trees_rs_rm$code[trees_rs$class==81] <- 'quag'
trees_rs_rm$code[trees_rs$class==74] <- 'plra'
trees_rs_rm$code[trees_rs$class==73] <- 'piun'
trees_rs_rm$code[trees_rs$class==66] <- 'pica'
trees_rs_rm$code[trees_rs$class==34] <- 'eugl'

ggplot() +
  geom_boxplot(data=dfyr2, aes(x="model", y=lai, fill='rhessys'),position=position_dodge2(padding=0.5)) + facet_wrap('code') + 
  geom_boxplot(data=trees_rs_rm, aes(x="data", y=lai, fill='alonzo'), position=position_dodge2(padding = 0.5)) + facet_wrap('code') + 
  ggtitle("compare LAI")

## Warning: Removed 2143 rows containing non-finite values (stat_boxplot).

ggplot() +
  geom_boxplot(data=dfyr2, aes(x="model", y=plantc, fill='rhessys')) + facet_wrap('code') + 
  geom_boxplot(data=trees_rs_rm, aes(x="data", y=carbon, fill='alonzo')) + facet_wrap('code', scales='free') + 
  ggtitle("compare plant carbon (kgC/m^2)")

## Warning: Removed 2143 rows containing non-finite values (stat_boxplot).

# ggplot() +
#   geom_boxplot(data=dfyr2, aes(x="model", y=height, fill='rhessys')) + facet_wrap('code') +
#   geom_boxplot(data=trees_rs, aes(x="data", y=height, fill='alonzo')) + facet_wrap('code') +
#   ggtitle("compare height (m)")

# ggplot() +
#   geom_boxplot(data=dfyr2, aes(x="model", y=rootdepth, fill='rhessys')) + facet_wrap('code')

Plots of NDVI and model outputs

# load in David's data 
#vis <- read.csv("/Volumes/GoogleDrive/My Drive/plant_params/DM_data/tree_and_turfgrass_monroe2_selected_weightedmeanVIs.csv")
vis <- read.csv("../../data/tree_and_turfgrass_monroe2_selected_weightedmeanVIs.csv")

trees_vi <- vis %>% dplyr::filter(sp_class==81|
                                    sp_class==74|
                                    sp_class==73|
                                    sp_class==66|
                                    sp_class==34)
vi_dens <- ggplot(trees_vi) + geom_density(aes(x=NDVI, col=as.factor(year))) + facet_wrap("sp_code")

rh_dens <- dfyr2 %>% ungroup() %>%
  dplyr::filter(year ==2011 | year==2014 | year==2017) %>% 
  ggplot() + geom_density(aes(x=lai,  col=as.factor(year))) + facet_wrap("code")

grid.arrange(vi_dens, rh_dens)

model_pc <- dfyr2 %>% ungroup() %>%
  dplyr::filter(year ==2011 | year==2014 | year==2017) %>% 
  pivot_wider(values_from=lai, names_from=year,  id_cols=c(run,code)) %>%
  mutate(decline = (`2014`-`2011`)/`2011`) %>% 
  mutate(recover = (`2017`-`2011`)/`2011`) %>%
  mutate(data="rhessys LAI") %>%
  dplyr::select(ID = run, sp_code=code, decline, recover, data)
model_pc$sp_code = toupper(model_pc$sp_code)

index_pc <- trees_vi %>% 
  dplyr::filter(year ==2011 | year==2014 | year==2017) %>% 
  pivot_wider(values_from=NDVI, names_from=year,  id_cols=c(ID,sp_code)) %>%
  mutate(decline = (`2014`-`2011`)/`2011`) %>% 
  mutate(recover = (`2017`-`2011`)/`2011`) %>% 
  mutate(data="NDVI") %>% 
  dplyr::select(-`2011`, -`2014`, -`2017`)

perc_changes <- rbind(model_pc, index_pc)

dec <- ggplot(perc_changes) + 
  geom_boxplot(aes(x=sp_code, y=decline, fill=data))

rec <- ggplot(perc_changes) + 
  geom_boxplot(aes(x=sp_code, y=recover, fill=data))

perc_changes %>% 
  pivot_longer(cols=c(decline, recover), names_to="change", values_to="perc_diff") %>% 
  ggplot() + 
  geom_hline(aes(yintercept=0), linetype="dashed") + 
  geom_boxplot(aes(x=change, y=perc_diff, col=data), fill='grey') + 
  facet_wrap('sp_code') + 
  theme_bw() + theme(legend.position="top") + 
  labs(y="percent difference to 2011 values") 

Resilience indices

Llorett et al. 2011 defines resilience as the capacity to reach pre-episode growth levels and relative resilience as the resilience weighted by the damage incurred during the episode.

By looking at time series as each year being relative the the pre-drought conditions, we can see how many years it takes to recover (or not at all). can also figure out which year to consider the drought occurring - where do we assign drought vs. post drought?

\(resistance = Drought / PreDrought\)

\(recovery = PostDrought / Drought\)

\(resilience = PostDrought / PreDrought\)

\(RelativeResilience = (PostDrought-Drought)/PreDrought\)

• this example uses coast live oak
• 2011 is year used as pre-drought

In the model, tree LAI responds the worst to the drought in 2015, after the months of experiencing the extreme heat in 2014 and 2015.

Plots of parameter distributions

p_dec <- read.csv("decid_params.csv")
p_ev <- readRDS("evg_broadleaf_params.rds")
p_con <- readRDS("pica_conifer_params.rds")

p_dec2 <- cbind(p_dec[,15:54], p_dec$code)
p_ev2 <- cbind(p_ev[,14:53], p_ev$code)
p_con2 <- cbind(p_con[,14:53], p_con$code)

p_names <- c(
      "pore_size_index",
      "psi_air",
      "run",
       "epc.leaf_cn",
       "epc.branch_turnover",
       "epc.gl_smax",
       "epc.flnr_age_mult",
       "epc.litter_moist_coef",
       "epc.psi_close",
       "mrc.q10",
       "epc.vpd_close",
       "epc.leaf_turnover",
       "epc.froot_turnover",
       "epc.storage_transfer_prop",
       "epc.height_to_stem_coef",
       "epc.waring_pa",
        "epc.root_distrib_parm",
       "epc.proj_sla",
       "epc.alloc_stemc_leafc",
       "epc.ext_coef",
       "specific_rain_capacity",
       "epc.flnr_sunlit",
       "epc.flnr_shade",
       "epc.netpabs_sunlit",
       "epc.netpabs_shade",
       "epc.netpabs_age_mult",
      "epc.cpool_mort_fract",
       "epc.min_percent_leafg",
      "epc.livewood_turnover",
      "epc.waring_pb",
      "epc.max_storage_percent",
      "epc.min_leaf_carbon",
      "epc.resprout_leaf_carbon",
      "epc.alloc_frootc_leafc",
      "epc.frootc_crootc",
      "epc.alloc_prop_day_growth",
      "epc.day_leafon",
       "epc.day_leafoff",
       "epc.ndays_expand",
       "epc.ndays_litfall", "code")

colnames(p_dec2) <- p_names
colnames(p_ev2) <- p_names
colnames(p_con2) <- p_names

all_p <- rbind(p_dec2, p_ev2, p_con2)

###### now can use to plot 

ggplot(all_p) + geom_boxplot(aes(x=as.factor(code), y=epc.proj_sla, fill=as.factor(code))) + ggtitle("proj_sla")

ggplot(all_p) + geom_boxplot(aes(x=as.factor(code), y=epc.leaf_turnover, fill=as.factor(code))) + ggtitle("leaf_turnover")

ggplot(all_p) + geom_boxplot(aes(x=as.factor(code), y=epc.waring_pa, fill=as.factor(code))) + ggtitle("waring_pa")

ggplot(all_p) + geom_boxplot(aes(x=as.factor(code), y=epc.waring_pb, fill=as.factor(code))) + ggtitle("waring_pb")

ggplot(all_p) + geom_boxplot(aes(x=as.factor(code), y=epc.branch_turnover, fill=as.factor(code))) + ggtitle("branch_turnover")

ggplot(all_p) + geom_boxplot(aes(x=as.factor(code), y=epc.froot_turnover, fill=as.factor(code))) + ggtitle("froot_turnover")

ggplot(all_p) + geom_boxplot(aes(x=as.factor(code), y=epc.root_distrib_parm, fill=as.factor(code))) + ggtitle("root distrib")

comb <- left_join(dfyr2, all_p, by=c('code','run')) %>%
  dplyr::filter(year==2011 | year==2014 | year==2017) %>% 
  group_by(year, run, code) %>% summarize_all(funs(mean), na.rm=T)

## Warning: Column `code` joining character vector and factor, coercing into
## character vector

key parameters for lai and plantc:

• needleleaf: pore size entry, psi air index, psi_close
• decid: waring_pa, proj_sla
• both: netpabs_sunlit, netpabs_shade, branch_turnover

ggplot(comb) + geom_jitter(aes(x=epc.branch_turnover, y=plantc, col=as.factor(year))) + facet_wrap('code') +
  ggtitle("branch_turnover")

## Warning: Removed 288 rows containing missing values (geom_point).

ggplot(comb) + geom_jitter(aes(x=epc.leaf_turnover, y=lai, col=as.factor(year))) + facet_wrap('code') +
  ggtitle("leaf_turnover") + geom_smooth(aes(x=epc.leaf_turnover, y=lai, col=as.factor(year)))

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## Warning: Removed 288 rows containing non-finite values (stat_smooth).

## Warning: Removed 288 rows containing missing values (geom_point).

ggplot(comb) + geom_jitter(aes(x=epc.waring_pa, y=lai, col=as.factor(year))) + facet_wrap('code') + ggtitle("waring pa")

## Warning: Removed 288 rows containing missing values (geom_point).

ggplot(comb) + geom_jitter(aes(x=epc.proj_sla, y=lai, col=as.factor(year))) + facet_wrap('code') + ggtitle("proj_sla")

## Warning: Removed 288 rows containing missing values (geom_point).

ggplot(comb) + geom_jitter(aes(x=epc.psi_close, y=lai, col=as.factor(year))) + facet_wrap('code') +
  ggtitle("psi_close")

## Warning: Removed 288 rows containing missing values (geom_point).